# installation 
Sys.setenv(http_proxy='https://10.118.113.5:3128')
Sys.setenv(https_proxy='https://10.118.113.5:3128')
Sys.setenv(ftp_proxy='https://10.118.113.5:3128')
Sys.setenv(http_proxy_user=ask)
Sys.setenv(https_proxy_user=ask)
Sys.setenv(ftp_proxy_user=ask)
source("http://www.bioconductor.org/biocLite.R")
biocLite("edgeR")


# start 
biocLite("edgeR")
library(edgeR)
# From the server

# part 1 read in data
setwd("/Users/naru/Box\ Documents/R\ code/data")   # error 1 
getwd()
raw.data <- read.table( file = "NC1.txt" , header = TRUE ) 
head( raw.data )

# part 2 
counts <- raw.data[ , -c(1,ncol(raw.data)) ]
rownames( counts ) <- raw.data[ , 1 ] # gene names
colnames( counts ) <- paste(c(rep("C_R",4),rep("T_R",3)),c(1:4,1:3),sep="") # sample names
head( counts )

# part 3 

dim( counts )
colSums( counts ) # Library Sizes
colSums( counts ) / 1e06 # Library Sizes in millions of reads
table( rowSums( counts ) )[ 1:30 ] # Number of genes with low counts
group <- c(rep("C", 4) , rep("T", 3))
cds <- DGEList( counts , group = group )
names( cds )
head(cds$counts) # original count matrix
cds$samples # contains a summary of your samples
sum( cds$all.zeros ) # How many genes have 0 counts across all samples
cds # or type the name of the object
cds <- cds[rowSums(1e+06 * cds$counts/expandAsMatrix(cds$samples$lib.size, dim(cds)) > 1) >= 3, ]
dim( cds )
cds <- calcNormFactors( cds )
cds$samples
# effective library sizes
cds$samples$lib.size * cds$samples$norm.factors
# To view the plot immediately
plotMDS.dge( cds , main = "MDS Plot for Count Data", labels = colnames( cds$counts ) )
# Output plot as a pdf
pdf( "MDS_plot_1_ex1.pdf" , width = 7 , height = 7 ) # in inches
plotMDS.dge( cds , main = "MDS Plot for Count Data", labels = colnames( cds$counts ) )
dev.off() # this tells [R] to close and stop writing to the pdf.

cds <- estimateCommonDisp( cds )
names( cds )
# The estimate
cds$common.dispersion
sqrt( 200 ) # poisson sd
sqrt( 200 + 200^2 * cds$common.dispersion ) # negative binomial sd
sqrt( 200 + 200^2 * cds$common.dispersion ) / sqrt( 200 ) # NB sd is over 2 times larger

# Default Setting
cds <- estimateTagwiseDisp( cds , prior.n = 10 )
names( cds )
summary( cds$tagwise.dispersion )
# More shrinkage/sqeezing toward the common
cds <- estimateTagwiseDisp( cds , prior.n = 25 )
summary( cds$tagwise.dispersion ) # not much changed, but the ends got squeezed in quite a bit.
# The recommended setting for this data set is the default of 10. Let’s stick with that.
cds <- estimateTagwiseDisp( cds , prior.n = 10 )
meanVarPlot <- plotMeanVar( cds , show.raw.vars=TRUE ,
                            show.tagwise.vars=TRUE ,
                            show.binned.common.disp.vars=FALSE ,
                            show.ave.raw.vars=FALSE ,
                            dispersion.method = "qcml" , NBline = TRUE ,
                            nbins = 100 ,
                            pch = 16 ,
                            xlab ="Mean Expression (Log10 Scale)" ,
                            ylab = "Variance (Log10 Scale)" ,
                            main = "Mean-Variance Plot" )


de.cmn <- exactTest( cds , common.disp = TRUE , pair = c( "C" , "T" ) )
de.tgw <- exactTest( cds , common.disp = FALSE , pair = c( "C" , "T" ) )
de.poi <- exactTest( cds , dispersion = 1e-06 , pair = c( "C" , "T" ) )
names( de.tgw )
de.tgw$comparison # which groups have been compared
head( de.tgw$table ) # results table in order of your count matrix.
head( cds$counts )


# Top tags for tagwise analysis
options( digits = 3 ) # print only 3 digits
topTags( de.tgw , n = 20 , sort.by = "p.value" ) # top 20 DE genes
# Back to count matrix for tagwise analysis
cds$counts[ rownames( topTags( de.tgw , n = 15 )$table ) , ]
# Sort tagwise results by Fold-Change instead of p-value
resultsByFC.tgw <- topTags( de.tgw , n = nrow( de.tgw$table ) , sort.by = "logFC" )$table
head( resultsByFC.tgw )
# Store full topTags results table
resultsTbl.cmn <- topTags( de.cmn , n = nrow( de.cmn$table ) )$table
resultsTbl.tgw <- topTags( de.tgw , n = nrow( de.tgw$table ) )$table
resultsTbl.poi <- topTags( de.poi , n = nrow( de.poi$table ) )$table
head( resultsTbl.tgw )

# Names/IDs of DE genes
de.genes.cmn <- rownames( resultsTbl.cmn )[ resultsTbl.cmn$adj.P.Val <= 0.05 ]
de.genes.tgw <- rownames( resultsTbl.tgw )[ resultsTbl.tgw$adj.P.Val <= 0.05 ]
de.genes.poi <- rownames( resultsTbl.poi )[ resultsTbl.poi$adj.P.Val <= 0.05 ]
# Amount significant
length( de.genes.cmn )
length( de.genes.tgw )
length( de.genes.poi )
# Percentage of total genes
length( de.genes.cmn ) / nrow( resultsTbl.cmn ) * 100
length( de.genes.tgw ) / nrow( resultsTbl.tgw ) * 100
length( de.genes.poi ) / nrow( resultsTbl.poi ) * 100
# Up/Down regulated summary for tagwise results
summary( decideTestsDGE( de.tgw , p.value = 0.05 ) ) # the adjusted p-values are used here


sum( de.genes.tgw %in% de.genes.cmn ) / length( de.genes.tgw ) * 100 # Tagwise to Common
sum( de.genes.cmn %in% de.genes.tgw ) / length( de.genes.cmn ) * 100 # Common to Tagwise
sum( de.genes.tgw %in% de.genes.poi ) / length( de.genes.tgw ) * 100 # Tagwise to Poisson
# Percent shared out of top 10, 100 & 1000 between tagwise and common
sum( de.genes.tgw[1:10] %in% de.genes.cmn[1:10] ) / 10 * 100
sum( de.genes.tgw[1:100] %in% de.genes.cmn[1:100] )
sum( de.genes.tgw[1:1000] %in% de.genes.cmn[1:1000] ) / 1000 * 100
# Percent shared out of top 10, 100 & 1000 between tagwise and poisson
sum( de.genes.tgw[1:10] %in% de.genes.poi[1:10] ) / 10 * 100
sum( de.genes.tgw[1:100] %in% de.genes.poi[1:100] )
sum( de.genes.tgw[1:1000] %in% de.genes.poi[1:1000] ) / 1000 * 100

par( mfrow=c(3 ,1) )
hist( resultsTbl.poi[de.genes.poi[1:100],"logConc"] , breaks=10 , xlab="Log Concentration" ,
      col="red" , xlim=c(-18,-6) , ylim=c(0,0.4) , freq=FALSE , main="Poisson: Top 100" )
hist( resultsTbl.cmn[de.genes.cmn[1:100],"logConc"] , breaks=25 , xlab="Log Concentration" ,
      col="green" , xlim=c(-18,-6) , ylim=c(0,0.4) , freq=FALSE , main="Common: Top 100" )
hist( resultsTbl.tgw[de.genes.tgw[1:100],"logConc"] , breaks=25 , xlab="Log Concentration" ,
      col="blue" , xlim=c(-18,-6) , ylim=c(0,0.4) , freq=FALSE , main="Tagwise: Top 100" )
par( mfrow=c(1,1) )

par( mfrow=c(2,1) )
plotSmear( cds , de.tags=de.genes.poi , main="Poisson" ,
           pair = c("C","T") ,
           cex = .35 ,
           xlab="Log Concentration" , ylab="Log Fold-Change" )
abline( h = c(-2, 2) , col = "dodgerblue" )
plotSmear( cds , de.tags=de.genes.tgw , main="Tagwise" ,
           pair = c("C","T") ,
           cex = .35 ,
           xlab="Log Concentration" , ylab="Log Fold-Change" )
abline( h=c(-2,2) , col="dodgerblue" )
par( mfrow=c(1,1) )


par( mfrow=c(2,1) )
plotSmear( cds , de.tags=de.genes.poi , main="Poisson" ,
           pair = c("C","T") ,
           cex = .35 ,
           xlab="Log Concentration" , ylab="Log Fold-Change" )
abline( h = c(-2, 2) , col = "dodgerblue" )
plotSmear( cds , de.tags=de.genes.tgw , main="Tagwise" ,
           pair = c("C","T") ,
           cex = .35 ,
           xlab="Log Concentration" , ylab="Log Fold-Change" )
abline( h=c(-2,2) , col="dodgerblue" )
par( mfrow=c(1,1) )


# Change column names to be specific to the analysis, logConc and logFC are the same in both.
colnames( resultsTbl.cmn ) <- c( "logConc" , "logFC" , "pVal.Cmn" , "adj.pVal.Cmn" )
colnames( resultsTbl.tgw ) <- c( "logConc" , "logFC" , "pVal.Tgw" , "adj.pVal.Tgw" )
# Below provides the info to re-order the count matrix to be in line with the order of the results.
wh.rows.tgw <- match( rownames( resultsTbl.tgw ) , rownames( cds$counts ) )
wh.rows.cmn <- match( rownames( resultsTbl.cmn ) , rownames( cds$counts ) )
head( wh.rows.tgw )
# Tagwise Results
combResults.tgw <- cbind( resultsTbl.tgw ,
                          "Tgw.Disp" = cds$tagwise.dispersion[ wh.rows.tgw ] ,
                          "UpDown.Tgw" = decideTestsDGE( de.tgw , p.value = 0.05 )[ wh.rows.tgw ] ,
                          cds$counts[ wh.rows.tgw , ] )
head( combResults.tgw )
# Common Results
combResults.cmn <- cbind( resultsTbl.cmn ,
                          "Cmn.Disp" = cds$common.dispersion ,
                          "UpDown.Cmn" = decideTestsDGE( de.cmn , p.value = 0.05 )[ wh.rows.cmn ] ,
                          cds$counts[ wh.rows.cmn , ] )
head( combResults.cmn )

wh.rows <- match( rownames( combResults.cmn ) , rownames( combResults.tgw ) )
combResults.all <- cbind( combResults.cmn[,1:4] ,
                          combResults.tgw[wh.rows,3:4] ,
                          "Cmn.Disp" = combResults.cmn[,5],
                          "Tgw.Disp" = combResults.tgw[wh.rows,5],
                          "UpDown.Cmn" = combResults.cmn[,6],
                          "UpDown.Tgw" = combResults.tgw[wh.rows,6],
                          combResults.cmn[,7:ncol(combResults.cmn)] )
head( combResults.all )
# Ouput csv tables of results
write.table( combResults.tgw , file = "combResults_tgw_ex1.csv" , sep = "," , row.names = TRUE )
write.table( combResults.cmn , file = "combResults_cmn_ex1.csv" , sep = "," , row.names = TRUE )
write.table( combResults.all , file = "combResults_all_ex1.csv" , sep = "," , row.names = TRUE )
